OVERVIEW

Credits I-Hsuan Lin, National Yang-Ming University Image credits: I-Hsuan Lin, National Yang-Ming University

OUR WORKFLOW

PASSAGGIO SOFTWARE FORMATO FILE INPUT FORMATO FILE OUTPUT FILE SUPPLEMENTARI
Controllo qualità delle reads FastQC .fastq (.fq; compresso: .fq.gz) .fastq (.fq; compresso: .fq.gz)
Rimozione adattatori e pulizia basata su qualità Trimmomatic .fastq (.fq; compresso: .fq.gz) .fastq (.fq; compresso: .fq.gz) File fasta con sequenze adattatori
Mappaggio delle reads STAR .fastq (.fq; compresso: .fq.gz) .sam / .bam (compresso) File fasta del genoma e file GTF dell’annotazione
Visualizzazione reads mappate IGV (Integrative Genomics Viewer) .bam (indicizzato .bai) Nessun output File del genoma e annotazione (già disponibile su IGV)
Controllo qualità allineamento (opzionale) Qualimap .bam Output tabulare
Conta delle reads featureCounts .sam / .bam .txt (testo) File GTF dell’annotazione
Analisi geni differenzialmente espressi DESeq2 / edgeR (pacchetti R, Bioconductor) .txt (testo) .txt / .csv (qualsiasi formato tabulare) Tabella con disegno sperimentale
Data mining, Gene Ontology Enrichment Analysis (GOEA) biomaRt, clusterprofiler oggetti R derivati da analisi DE grafici

DAY 1 - 22 maggio 2018

RAW READS MANIPULATION

Accendere il computer e aprire il terminale (ctrl + alt + T)

Spostarsi nella cartella Scaricati, create una nuova cartella chiamandola tc

cd Scaricati/
mkdir tc

Entrare nella cartella appena creata e valutarne il contenuto della cartella digitare il seguente comando

cd tc/
ls -lrth

La cartella dovrebbe essere vuota. Come prima cosa dobbiamo recuperare i file grezzi (FASTQ) del nostro esperimento. Per recuperare i file ci connettiamo al nostro server di laboratorio e con il protocollo di trasferimento ssh copiamo in locale (la cartella dove siamo in questo momento) i 2 file corrispondenti alle reads forward e alle reads reverse. Lavoreremo con un solo campione, ma come sapete l’esperimento consiste di 3 repliche biologiche per i non trattati e 3 repliche biologiche per quelli trattati. I comandi che qui eseguiremo possono essere ripetuti su più file contemporaneamente tramite appositi comandi come i cicli for

# prima scarichiamo il primo file (attenzione al punto in fondo al comando, indica current directory)
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/Sample_Untreated_1_1_subset.fastq.gz .
# inserire la password quando ce la chiede ()
# controlliamo che sia stato scaricato
ls -lrth
# ora scarichiamo il secondo file
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/Sample_Untreated_1_2_subset.fastq.gz .
# inserire nuovamente la password
# controlliamo
ls -lrth

Come potete osservare la dimensione di questi file è molto bassa rispetto alla dimensione tipica di un file che risulta da un esperimento di sequenziamento come lo abbiamo descritto. Questo perché i file originali sono stati copiati e ridotti in dimensioni con lo scopo di riuscire a manipolarli e lavorarci nei tempi della lezione, principalmente per il fatto che le capacità dei computer che stiamo usando potrebbero essere limitate e quindi non in grado di eleaborare questi file in tempi brevi

Iniziamo con esplorare questi file

# "apriamo" il file
less -S Sample_Untreated_1_1_subset.fastq.gz
# contiamo il numero di reads
zcat Sample_Untreated_1_1_subset.fastq.gz | echo $((`wc -l`/4))
zcat Sample_Untreated_1_2_subset.fastq.gz | echo $((`wc -l`/4))
# il numero di reads F e R coincide?
FASTQ file

FASTQ file

QUALITY CONTROL

Per effettuare il controllo di qualità sulle nostre reads utilizzeremo il programma FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Questo software (scritto in Java) prende come input un numero a piacimento di reads e calcola varie statistiche che possono dare un’idea della qualità dei nostri dati

# creiamo una nuova cartella dove scaricheremo il programma e salveremo i risultati
mkdir fastqc
cd fastqc/
# scarichiamo il programma con il comando wget
wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip
# scompattiamo la cartella
unzip fastqc_v0.11.7.zip
# rimuoviamo cartella compressa
rm fastqc_v0.11.7.zip
# vediamo cosa c'è dentro la cartella senza spostarci dentro (percorso relativo alla nostra posizione)
ls -lrth FastQC/
# il file che ci interessa, cioè il programma che andremo a richiamare è "fastqc"
# dobbiamo prima renderlo eseguibile
chmod 755 FastQC/fastqc
# controlliamo che funzioni
./FastQC/fastqc -h

A questo punto possiamo lanciare il programma dando come input i due file con le reads F e R. Il comando è molto semplice, basta specificare i nomi dei file che vogliamo controllare, e la cartella dove vogliamo che vengano salvati i report coi risultati (la cartella deve già esistere quindi prima la creiamo)

# creiamo cartella per salvare i risultati
mkdir results
# ci spostiamo una cartella su, dove stanno i nostri file
cd ..
# eseguiamo il programma (dobbiamo specificare il percorso dell'eseguibile rispetto a dove siamo) e specifichiamo inoltre il percorso della cartella di output
./fastqc/FastQC/fastqc Sample_Untreated_1_1_subset.fastq.gz Sample_Untreated_1_2_subset.fastq.gz -o fastqc/results/

Per ogni file di input vengono creati 2 file di output, un file html e una cartella compressa. Andiamo ad esplorare il report in formato html aprendolo con Firefox

FastQC report

FastQC report

A questo punto commentiamo insieme i risultati del controllo qualità confrontandoci con qualità di riferimento sul sito di FastQC e con la spiegazione dei moduli di analisi (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/)

Phred quality

Phred quality

Dopo aver deciso insieme quali parametri utilizzare per correggere la qualità delle nostre reads ed aver identificato gli adattatori da rimuovere, andiamo ad effettuare la rimozione degli adattatori e il quality trimming

ADAPTERS REMOVAL AND QUALITY TRIMMING

Per effettuare queste due operazioni useremo un unico programma che si chiama Trimmomatic (http://www.usadellab.org/cms/?page=trimmomatic). Il suo utilizzo è molto simile a quello di FastQC (anche Trimmomatic è scritto in Java) e necessita come input le sequenze grezze. Per ogni file di input (nel nostro caso 2 file) verranno creati 2 file di output, uno indicato come paired, l’altro come unpaired; è necessario indicare il nome di ciascuno dei file di output. Le opzioni di Trimmomatic sono:

  • PE indica che lavoriamo su file paired-end
  • -trimlog indica che vogliamo salvare un file di log
  • -threads indica il numero di processori che vogliamo utilizzare per lanciare il programma
  • ILLUMINACLIP:ALL_ADAPTERS.fa:2:30:10 è la parte che indica la rimozione degli adattatori. Il file ALL_ADAPTERS.fa è il file in formato fasta che contiene la sequenza degli adattatori maggiormente utilizzati
  • LEADING:16 indica la rimozione di N basi in testa alla sequenza, sotto la qualità indicata
  • TRAILING:16 indica la rimozione di N basi in coda alla sequenza, sotto la qualità indicata
  • SLIDINGWINDOW:4:25 effettua un trimming basato su una “finestra” di basi, rimuovendo le basi nel momento in cui la qualità media per quella finestra scende sotto la qualità indicata
  • MINLEN:36 indica la rimozione delle reads che sono più corte della lunghezza indicata
# ci spostiamo nella cartella tc, creiamo nuova cartella e scarichiamo il programma
# usiamo percorso assoluto da qualsiasi punto del computer per spostarci dentro una qualsiasi cartella
cd /home/usr/Scaricati/tc # usr è specifico per il vostro utente (il mio è dlfptr70)
mkdir Trimmomatic
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.36.zip -O Trimmomatic/Trimmomatic-0.36.zip # con l'opzione -O specifichiamo dove vogliamo scaricare il file e con che nome
# decomprimiamo senza cambiare cartella e rimuoviamo l'archivio
unzip Trimmomatic/Trimmomatic-0.36.zip -d Trimmomatic/ && rm Trimmomatic/Trimmomatic-0.36.zip
# controlliamo il funzionamento del programma
java -jar Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar -h # -h è l'opzione per help

A questo punto possiamo lanciare Trimmomatic sui nostri file; anche se dal controllo qualità non abbiamo identificato sequenze di adattatori come contaminanti, un classico approccio potrebbe essere quello di dare come input di file degli adattatori il file fornito dal programma che quindi rimuoverà tutte le sequenze corrispondenti a quelle che trova

# creiamo file unendo in uno tutti gli adattatori
cat Trimmomatic/Trimmomatic-0.36/adapters/*.fa > Trimmomatic/Trimmomatic-0.36/adapters/ALL_ADAPTERS.fa
# cat è un comando che unisce file "row-wise"; l'asterisco (*.fa) indica di selezionare tutti i file che hanno estensione .fa
java -jar Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 10 -trimlog \
Trimmomatic/log_trimmomatic \
Sample_Untreated_1_1_subset.fastq.gz Sample_Untreated_1_2_subset.fastq.gz \
Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz Trimmomatic/Sample_Untreated_1_1_subset_trim_unpaired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz Trimmomatic/Sample_Untreated_1_2_subset_trim_unpaired.fastq.gz \
ILLUMINACLIP:Trimmomatic/Trimmomatic-0.36/adapters/ALL_ADAPTERS.fa:2:30:10 \
LEADING:16 TRAILING:16 SLIDINGWINDOW:4:25 MINLEN:36 &> Trimmomatic/log_trimmomatic_screen
# con il simbolo &> andiamo a salvare su un file anche l'output che normalmente sarebbe stampato sullo schermo

Dovrebbe terminare dopo qualche secondo senza aver stampato a schermo alcuna informazione. Andiamo a vedere se ha creato i file ed esaminiamo i file di log (log_trimmomatic e log_trimmomatic_screen); il primo contiene le seguenti info per ciascuna read:

  • the read name
  • the surviving sequence length
  • the location of the first surviving base, aka. the amount trimmed from the start
  • the location of the last surviving base in the original read
  • the amount trimmed from the end

Il secondo contiene le statistiche generali

# guardiamo prima le statistiche generali
less -S Trimmomatic/log_trimmomatic_screen
# poi l'altro file
less -S Trimmomatic/log_trimmomatic

A questo punto ripetiamo il controllo di qualità, questa volta sulle reads trimmate, solo sui file paired

# creiamo cartella dove salveremo nuovo controllo qualità
mkdir fastqc/results_after_trimming
# lanciamo FastQC sulle sequenze trimmate
./fastqc/FastQC/fastqc Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz -o fastqc/results_after_trimming/

Andiamo ad esplorare i report e vedere se le statistiche sulla qualità sono cambiate

MAPPING

Ora i nostri file sono pronti per essere allineati sul genoma di riferimento (Arabidopsis). Per quasi tutti i programmi che effettuano allineamento il primo passaggio consiste della creazione dell’indice del genoma. Perché viene creato l’indice del genoma: citazione presa da qualcuno ma non ricordo chi (modificata) “Se entro in una biblioteca con la frase “La genetica è il ramo della biologia che si occupa del materiale ereditario” e voglio trovare il libro in cui è scritta, farò molta fatica a trovarlo, a meno che all’ingresso della biblioteca non si trovi un catalogo che mi indirizza almeno verso la sezione “Scienza””. L’indice quindi consiste in una serie di file che riassumono le principali caratteristiche del genoma di un organismo. Per effettuare l’allineamento useremo un programma chiamato STAR (https://github.com/alexdobin/STAR). Andiamo a scaricarlo insieme ai file aggiornati del genoma e del trascrittoma di Arabidopsis

# siamo nella cartella tc
mkdir STAR
# scarichiamo e salviamo nella cartella appena creata
wget https://github.com/alexdobin/STAR/archive/2.6.0c.tar.gz -O STAR/2.6.0c.tar.gz
# decomprimiamo e rimuoviamo archivio
tar -xf STAR/2.6.0c.tar.gz -C STAR/ && rm STAR/2.6.0c.tar.gz
# controlliamo il funzionamento
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR -h

Per creare l’indice del genoma, abbiamo bisogno del genoma. Andiamo a scaricarlo dal sito di Ensembl Plants

# genoma
mkdir genome
wget ftp://ftp.ensemblgenomes.org/pub/release-39/plants/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz -O genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
# trascrittoma
wget ftp://ftp.ensemblgenomes.org/pub/release-39/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.39.gtf.gz -O genome/Arabidopsis_thaliana.TAIR10.39.gtf.gz
# i file devono essere decompressi
gunzip genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gunzip genome/Arabidopsis_thaliana.TAIR10.39.gtf.gz

A questo punto possiamo creare l’indice con il comando di STAR --runMode genomeGenerate. Vediamo nel dettaglio le altre opzioni:

  • ./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR chiamiamo il programma
  • --runThreadN numero di threads con cui vogliamo lanciare il programma
  • --runMode indica che vogliamo creare l’indice
  • --genomeDir indica la cartella dove vogliamo che venga salvato l’output
  • --genomeFastaFiles cartella dove è presente il file del genoma
  • --sjdbGTFfile indichiamo che vogliamo usare anche il trascrittoma per aiutarlo nella creazione dell’indice
  • --sjdbOverhang specifica la lunghezza delle reads meno 1. Serve per creare il database delle giunzioni di splicing

Lanciamo il comando con le opzioni e i percorsi adattati al nostro caso

# creiamo cartella per l'indice all'interno della cartella del genoma
mkdir -p genome/index
# lanciamo la creazione dell'indice
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR \
--runThreadN 10 \
--runMode genomeGenerate \
--genomeDir genome/index \
--genomeFastaFiles genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa \
--sjdbGTFfile genome/Arabidopsis_thaliana.TAIR10.39.gtf \
--sjdbOverhang 100

Con 10 core dovrebbe impiegare circa 3 minuti. Andiamo a vedere se ha creato i file nella cartella index

# spostiamo il file Log.out dove dovrebbe stare
mv Log.out genome/index/
# guardiamo i file appena creati
ls -lrht genome/index
# raccolgono diverse statistiche sul genoma di Arabidopsis, comprese info sul db di splicejunctions

Possiamo ora andare ad effettuare l’allineamento. Ovviamente per l’allineamento useremo le reads trimmate e pulite (le paired). Vediamo nel dettaglio i comandi:

  • ./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR chiamiamo il programma
  • --runThreadN numero di threads con cui vogliamo lanciare il programma
  • --genomeDir dove sta la cartella con l’indice
  • --readFilesIn qui scriviamo il percorso e il nome dei campioni da mappare
  • --readFilesCommand qual è il formato di comrpessione dei file da mappare
  • --outFileNamePrefix percorso e prefisso che vogliamo assegnare ai file di output
  • --alignIntronMax lunghezza massima (stimata) degli introni del nostro organismo
  • --outSAMtype tipo es estensione del file di output

Attenzione che la maggior parte dei programmi di allineamento sono creati con impostazioni di default basate sul genoma umano –> aggiustare i parametri sulla base del proprio organismo! Andiamo a lanciare il comando

# creiamo cartella dove salveremo file di allineamento
mkdir -p STAR/aligned
# lanciamo l'allineamento
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR \
--runThreadN 10 \
--genomeDir genome/index/ \
--readFilesIn Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz \
--readFilesCommand zcat \
--outFileNamePrefix STAR/aligned/SU_1_ \
--alignIntronMax 2000 \
--outSAMtype BAM SortedByCoordinate

Se è tutto a posto l’allineamento dovrebbe impiegare circa 10 secondi. Come vedete la lunghezza massima media per gli introni di Arabidopsis è di circa 2000 bp, e gli ho detto che voglio che il file di output sia in formato .bam (formato allineamento binario). Il file .bam è in un formato che non può essere aperto con il terminale, ma che andremo poi a visualizzare con un apposito software. Per completezza vi inserisco un immagine dello stesso file di allineamento ma in un formato visualizzabile e consultabile sul terminale. Il formato di questo file è .sam (Sequence Alignment Format), formato non molto utilizzato perché non essendo compresso occupa molto spazio

SAM file

SAM file

La prima riga descrive il numero della versione del file

Dalla seconda alla ottava sono riportate le lunghezze in paia di basi dei cromosomi di Arabidopsis

La nona e la decima riga indicano il programma utilizzato per l’allineamento e il comando digitato

Dalla undicesima riga ogni riga coincide con una read, le colonne sono descritte di seguito

La colonna più interessante è forse la sesta (CIGAR) dove vengono specificati nel dettaglio gli allineamenti di ciascuna reads con il reference

SAM file description

SAM file description

Andiamo ora a vedere i nostri file creati

# il file che ci interessa è chiamato SU_1_Aligned.sortedByCoord.out.bam
ls -rht STAR/aligned/
# leggiamo qualche statistica in questo altro file
less -S STAR/aligned/SU_1_Log.final.out

Quanto vi risulta la percentuale di Uniquely mapped reads?

ALIGNMENT VISUALIZATION

Dopo aver effettuato l’allineamento può essere utile effettuare un controllo visivo sull reads allineate. A tal fine si utilizza spesso un software chiamato IGV (http://software.broadinstitute.org/software/igv/) che permette tra le altre cose di visualizzare i file dell’allineamento dopo aver caricato un genoma di riferimento

Integrative Genomics Viewer

Integrative Genomics Viewer

Andiamo a scaricare il programma dal sito http://software.broadinstitute.org/software/igv/download cliccando in basso su “Launch with 2 GB”. Attendiamo lo scaricamento e apriamo il programma con doppio clic. Una volta aperto in alto a sinistra selezionare il genoma di A. thaliana TAIR10. A questo punto possiamo caricare i risultati del nostro allineamento, non prima però di aver creato un indice per il file .bam. Un altro indice direte voi, sì questa volta però non eseguiremo nessun passaggio ma vi farò prelevare il file dell’indice già preparato dal nostro server. Per completezza vi mostro il comando utilizzato per creare l’indice con Samtools

samtools index -b SU_1_Aligned.sortedByCoord.out.bam

Su questi computer non possiamo installare il programma che si utilizza, tra le altre cose, anche per la creazione degli indici per i file di allineamento. Il programma si chiama Samtools (http://www.htslib.org/) ed è molto utilizzato in bioinformatica perché svolge numerose funzioni. Andiamo a copiare nella cartella dove sta il file .bam appena creato da noi, il file dell’indice, che permette che il file .bam venga visualizzato correttamente in IGV

scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam.bai STAR/aligned/
# il file che vedete deve avere lo stesso nome del file di allineamento, ma ha un'estensione diversa, .bai

A questo punto potete tornare su IGV e cliccare su “File -> Load from file” e selezionare dalla cartella aligned il file SU_1_Aligned.sortedByCoord.out.bam. Selezioniamo Chr1 e zoomiamo sino a che la rotella inizia a girare, a quel punto possiamo visualizzare le nostre reads allineate ai geni del genoma di riferimento

Non sarà così facile trovare un po’ di reads!

Non sarà così facile trovare un po’ di reads!

Nel caso non funzionasse, copiare dal server alla vostra cartella, sia il file .bam che il .bai (potrebbe essere che il bai debba essere creato dallo stesso bam per funzionare la visualizzazione)

# bam
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam STAR/aligned/
# bai
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam.bai STAR/aligned/

La lunghezza degli introni è stata rispettata?

DAY 2 - 29 maggio 2018

READS COUNTING

Una volta controllato che l’allineamento è andato a buon fine e siamo soddisfatti di come il programma di mapping ha lavorato, possiamo spostarci nel prossimo step verso l’analisi dei geni, ovvero la conta delle reads. Per conta delle reads si intende il contare quante reads sono state assegnate (hanno “mappato”) a un determinato gene (o trascritto, o esone, le cosiddette genomic features) in seguito al mapping. Per questo passaggio utilizzeremo il software featureCounts (http://subread.sourceforge.net/) del pacchetto Subread. Scarichiamo e testiamo il programma

# nella cartella tc/ creiamo cartella
mkdir featureCounts
# scarichiamo dal server perché più comodo
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/featureCounts featureCounts/
# testiamo funzionamento
./featureCounts/featureCounts -v

Vediamo le opzioni e i parametri di featureCounts:

  • featureCounts lanciamo il programma
  • -a forniamo il file di annotazione del trascrittoma
  • -o percorso e nome file di output
  • -p indichiamo che lavoriamo su reads paired e che verranno contati i frammenti e non le reads
  • -t indica quale feature contare
  • -g indica che vogliamo raggrppare tutte le feature in questa meta-feature
  • -T quanti threads vogliamo usare per correre il programma
  • in.bam nome del file di allineamento che vogliamo come input

A questo punto andiamo a lanciare il programma, ricordandoci che nel nostro caso vogliamo contare le feature del file di allineamento .bam

# creiamo cartella per output
mkdir -p featureCounts/output
# lanciamo la conta
./featureCounts/featureCounts \
-a genome/Arabidopsis_thaliana.TAIR10.39.gtf \
-o featureCounts/output/SU_1_counts.txt \
-p \
-t exon \
-g gene_id \
-T 10 \
STAR/aligned/SU_1_Aligned.sortedByCoord.out.bam

Dovrebbe terminare in qualche secondo. Andiamo ad esplorare i file di output

ls -lrht featureCounts/output/
less -S featureCounts/output/SU_1_counts.txt.summary
less -S featureCounts/output/SU_1_counts.txt
# vediamo solo le info che ci interessano
cut -f 1,7 featureCounts/output/SU_1_counts.txt | sed 1d | less -S

Le colonne che ci interessano sono la prima e l’ultima (7). Perché per alcuni geni abbiamo più di un valore per le colonne 2, 3, 4 e 5? I passaggi che abbiamo visto finora in una situazione normale andrebbero ripetuti per tutti i campioni che abbiamo fatto sequenziare, sino ad ottenere una matrice di conte, come l’ultimo file che abbiamo visto, per ciascun campione. Da questa matrice di conte (counts matrix) creeremo un file unico che sarà l’input per la fase successiva. Visto che nella nostra esercitazione abbiamo elaborato solo un file, vi fornirò io la counts matrix completa (6 campioni, 3 controlli e 3 trattati)

# scarichiamo dal server il file con le conte per i 3 controlli e i 3 trattati
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/countDF_unvs3h featureCounts/
# vediamo il file
less -S featureCounts/countDF_unvs3h

INSTALLING R AND RSTUDIO

A questo punto abbiamo concluso la prima fase del corso che abbiamo svolto usando solamente il terminale. Per la seconda fase utilizzeremo invece il software R e il suo ambiente di sviluppo integrato (IDE) chiamato RStudio. Per prima cosa vi invito a scaricare e installare R ed RStudio seguendo le istruzioni di seguito a seconda del vostro SO:

Per qualsiasi problema scrivetemi all’indirizzo delfino.pietro@gmail.com o possiamo vedere domani l’installazione in aula, ma sarebbe meglio arrivare già con entrambi i programmi installati e alcuni pacchetti che elenco di seguito, con i rispettivi comandi per l’installazione. Aprire RStudio e nella console digitare i seguenti comandi con sfondo grigio, in ordine, uno per uno:

  • Bioconductor
    • source("https://bioconductor.org/biocLite.R")
    • biocLite()
  • tidyverse
    • install.packages("tidyverse")
  • DESeq2
    • source("https://bioconductor.org/biocLite.R")
    • biocLite("DESeq2")
  • edgeR
    • source("https://bioconductor.org/biocLite.R")
    • biocLite("edgeR")
  • biomaRt
    • source("https://bioconductor.org/biocLite.R")
    • biocLite("biomaRt")
  • clusterProfiler
    • source("https://bioconductor.org/biocLite.R")
    • biocLite("clusterProfiler")
  • pheatmap
    • install.packages("pheatmap")
  • RColorBrewer
    • install.packages("RColorBrewer")
  • Genome wide annotation of Arabidopsis
    • source("https://bioconductor.org/biocLite.R")
    • biocLite("org.At.tair.db")

DIFFERENTIAL GENE EXPRESSION ANALYSIS

DESEQ2

Torniamo alla lezione. Apriamo RStudio e ci spostiamo nella cartella dove abbiamo il file completo con le conte

setwd("/home/usr/Scaricati/tc/featureCounts") 
# "usr" è variabile a seconda del vostro nome utente, usate il tab per completare il percorso
# carichiamo la libreria
library(DESeq2)

# carichiamo il file
cts <- as.matrix(read.table("featureCounts/countDF_unvs3h", sep="\t",
                            header=T, col.names = c("ID","un1","un2","un3","NO3_1","NO3_2","NO3_3"),
                            row.names = "ID"))
# non eseguire questo comando, è solo per me
cts <- as.matrix(read.table("Scaricati/tc/featureCounts/countDF_unvs3h", sep="\t",
                           header=T, col.names = c("ID","un1","un2","un3","NO3_1","NO3_2","NO3_3"),
                           row.names = "ID"))

# vediamo cosa abbiamo caricato
head(cts)
##           un1  un2  un3 NO3_1 NO3_2 NO3_3
## AT1G01010 248  535  190   830   839   497
## AT1G01020 255  684  943   842   679   909
## AT1G01030 546  607  568    61    45    84
## AT1G01040 690 1970 1386   607   687   419
## AT1G01046   2    7    7     1     3     0
## AT1G01050 597 1252 2008  1015  1054  1061
# creiamo oggetto per il disegno sperimentale
coldata <- data.frame(row.names = colnames(cts))
coldata$condition <- c(rep("untreated", 3), rep("treated", 3))
coldata$time <- c(rep("0", 3), rep("3", 3))
str(coldata)
## 'data.frame':    6 obs. of  2 variables:
##  $ condition: chr  "untreated" "untreated" "untreated" "treated" ...
##  $ time     : chr  "0" "0" "0" "3" ...
# convertiamo le colonna in factor per poterle usare nella formula
coldata$condition <- as.factor(coldata$condition)
coldata$time <- as.factor(coldata$time)
str(coldata)
## 'data.frame':    6 obs. of  2 variables:
##  $ condition: Factor w/ 2 levels "treated","untreated": 2 2 2 1 1 1
##  $ time     : Factor w/ 2 levels "0","3": 1 1 1 2 2 2
# costruiamo l'oggetto DESeq
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ condition )

# ordiniamo i livelli per usare "untreated" come confronto
dds$condition <- relevel(dds$condition, ref = "untreated")

# vediamo info sull'oggetto creato
dds
## class: DESeqDataSet 
## dim: 33602 6 
## metadata(1): version
## assays(1): counts
## rownames(33602): AT1G01010 AT1G01020 ... ATMG01400 ATMG01410
## rowData names(0):
## colnames(6): un1 un2 ... NO3_2 NO3_3
## colData names(2): condition time
# rimuoviamo geni mai espressi
dds <- dds [rowSums(counts(dds)) > 0, ]

dds
## class: DESeqDataSet 
## dim: 24383 6 
## metadata(1): version
## assays(1): counts
## rownames(24383): AT1G01010 AT1G01020 ... ATMG01400 ATMG01410
## rowData names(0):
## colnames(6): un1 un2 ... NO3_2 NO3_3
## colData names(2): condition time
# quanti erano?

Data exploration with DESeq2

library(tidyverse)
# trasformiamo i dati per la  PCA; correggiamo per rendere dataset omoschedastico
# facciamo 2 trasformazioni proprie di DESeq2, rlog e vst
rld <- rlog(dds, blind = FALSE)
head(assay(rld),3)
##                un1      un2      un3    NO3_1    NO3_2    NO3_3
## AT1G01010 8.878895 8.722974 7.173810 9.572869 9.646944 8.948716
## AT1G01020 8.986913 9.101649 9.220015 9.659787 9.448712 9.779084
## AT1G01030 9.716037 8.717315 8.348288 6.198081 5.911402 6.591842
vsd <- vst(dds, blind = FALSE)
head(assay(vsd),3)
##                 un1      un2      un3    NO3_1    NO3_2    NO3_3
## AT1G01010  8.992480 8.825655 7.253611 9.743679 9.824669 9.067229
## AT1G01020  9.029901 9.155565 9.283881 9.763533 9.532699 9.894374
## AT1G01030 10.075139 8.994522 8.605132 6.542427 6.302154 6.885296
# calcoliamo PCA sui dati trasformati con rlog usando la funzione di R,
# dist, che calcola distanze Euclidee sulla matrice trasposa rld
sampleDists <- dist(t(assay(rld)))
sampleDists
##             un1       un2       un3     NO3_1     NO3_2
## un2    58.79275                                        
## un3    75.19080  64.52936                              
## NO3_1 199.18353 195.74731 200.74712                    
## NO3_2 197.94528 194.40729 198.29539  45.88588          
## NO3_3 185.94000 186.16615 177.45544  93.93457  93.37848
# plottiamo heatmap delle distanze specificando noi stessi la distanza che abbiamo
# appena calcolato con l'argomento "clustering distance"
sampleDistMatrix <- as.matrix( sampleDists )
rownames(sampleDistMatrix) <- rownames(coldata)
colnames(sampleDistMatrix) <- rownames(coldata)
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors
         )

# plottiamo PCA
# custom PCA plot
pcaData <- plotPCA(rld, intgroup=c("condition", "time"), returnData=TRUE)
pcaData.vsd <- plotPCA(vsd, intgroup=c("condition", "time"), returnData=TRUE)

ggplot(pcaData, aes(PC1, PC2, color=condition, shape=time)) +
  geom_text_repel(aes(label=pcaData$name)) +
  geom_point(size=3) +
  theme_bw() +  ggtitle("Data transformed with rld")

ggplot(pcaData.vsd, aes(PC1, PC2, color=condition, shape=time)) +
  geom_text_repel(aes(label=pcaData$name)) +
  geom_point(size=3) +
  theme_bw() +  ggtitle("Data transformed with vsd")

Dopo alcune analisi esplorative possiamo eseguire il test della differenza di espressione tra i geni vero e proprio. Teniamo presente che DESeq2 lavora sulle conte grezze per effettuare il test, e non sui dati trasformati, quindi useremo l’oggetto dds e non le trasformazioni rld o vsd

# comando che effettua l'analisi vera e propria
dds <- DESeq(dds)

# salviamo oggetto con risultati specificando quali confronti
res <- results(dds, contrast=c("condition","treated","untreated"), alpha = 0.01)
summary(res)
## 
## out of 24383 with nonzero total read count
## adjusted p-value < 0.01
## LFC > 0 (up)     : 4570, 19% 
## LFC < 0 (down)   : 4445, 18% 
## outliers [1]     : 11, 0.045% 
## low counts [2]   : 3782, 16% 
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# we can order our results table by the smallest p-value
resOrdered <- res[order(res$pvalue),]
resOrdered
## log2 fold change (MLE): condition treated vs untreated 
## Wald test p-value: condition treated vs untreated 
## DataFrame with 24383 rows and 6 columns
##             baseMean log2FoldChange     lfcSE       stat        pvalue
##            <numeric>      <numeric> <numeric>  <numeric>     <numeric>
## AT1G23740   5381.069      -5.994458 0.1754464  -34.16690 7.503470e-256
## AT1G22400  24333.995       7.492121 0.2250072   33.29724 4.232379e-243
## AT2G04040   9615.186       9.989993 0.3080733   32.42733 1.130964e-230
## AT2G21640   8584.189       6.425284 0.1986819   32.33956 1.945127e-229
## AT2G04050   8942.769       9.455977 0.3116858   30.33817 3.598939e-202
## ...              ...            ...       ...        ...           ...
## AT4G25470  773.59979      0.1412803  1.313684  0.1075451            NA
## AT5G35935   35.60476      4.6814514  2.768811  1.6907802            NA
## AT5G37770 1252.46550      0.2171257  1.240675  0.1750061            NA
## AT5G44430  872.74506     -1.0961080  2.117555 -0.5176291            NA
## AT5G62920  737.38478     -1.2402995  1.353632 -0.9162750            NA
##                    padj
##               <numeric>
## AT1G23740 1.544964e-251
## AT1G22400 4.357235e-239
## AT2G04040 7.762185e-227
## AT2G21640 1.001254e-225
## AT2G04050 1.482043e-198
## ...                 ...
## AT4G25470            NA
## AT5G35935            NA
## AT5G37770            NA
## AT5G44430            NA
## AT5G62920            NA

Volcano plot

# Threshold on the adjusted p-value
alpha <- 0.01
# colors for the points
cols <- densCols(res$log2FoldChange, -log10(res$pvalue))
# plot
plot(res$log2FoldChange, -log10(res$padj), col=cols, panel.first=grid(),
     main="Volcano plot", xlab="Effect size: log2(fold-change)", ylab="-log10(adjusted p-value)",
     pch=20, cex=0.6)
abline(v=0)
abline(v=c(-2,2), col="brown")
abline(h=-log10(alpha), col="brown")

#gn.selected <- abs(res$log2FoldChange) > 5 & res$padj < alpha 

#text(res$log2FoldChange[gn.selected],
#     -log10(res$padj)[gn.selected],
#     lab=rownames(res)[gn.selected ], cex=0.4)

Heatmap

# plottiamo con una heatmap i 20 geni con maggior varianza tra le condizioni
# sulla base della trasformazione rlog
topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 20)
mat  <- assay(rld)[ topVarGenes, ]
head(mat)
##                un1      un2      un3    NO3_1    NO3_2    NO3_3
## AT1G05680 7.944321 7.092012 7.958417 16.34366 16.43114 15.78029
## AT2G04040 5.576480 5.085690 5.407554 13.74236 13.96226 13.52739
## AT1G66570 4.798428 3.583184 4.159822 12.59125 12.24675 11.72936
## AT2G04050 5.902226 5.446585 5.604259 13.71389 13.92276 13.29248
## AT5G24640 5.522707 4.427213 4.520819 12.62856 12.79483 12.80844
## AT4G12735 6.341879 5.004217 6.125653 13.73173 13.67158 13.75279
mat  <- mat - rowMeans(mat)
head(mat)
##                 un1       un2       un3    NO3_1    NO3_2    NO3_3
## AT1G05680 -3.980652 -4.832960 -3.966555 4.418686 4.506166 3.855316
## AT2G04040 -3.973809 -4.464599 -4.142735 4.192072 4.411967 3.977103
## AT1G66570 -3.386370 -4.601614 -4.024977 4.406453 4.061948 3.544560
## AT2G04050 -3.744808 -4.200449 -4.042775 4.066861 4.275726 3.645446
## AT5G24640 -3.261054 -4.356549 -4.262943 3.844800 4.011069 4.024677
## AT4G12735 -3.429430 -4.767092 -3.645656 3.960424 3.900273 3.981482
anno <- as.data.frame(colData(rld)[, c("condition","time")])
pheatmap(mat, annotation_col = anno)

Scatter plot delle conte

# plottiamo le conte normalizzate del gene con il padj minore; o fold change minore
plotCounts(dds, gene=which.min(res$padj), intgroup="condition")

plotCounts(dds, gene=which.min(res$log2FoldChange), intgroup="condition")

# se siamo interessati a un gene in particolare possiamo scriverne il nome
plotCounts(dds, gene="AT1G66570", intgroup="condition")

GENE ONTOLOGY ANALYSIS

Gene ontology: Il discorso sull’essere dei geni

www.geneontology.org

www.geneontology.org

Per esplorare le gene ontology dei geni che abbiamo trovato come differenzialmente espressi utilizzeremo il pacchetto clusterprofiler (https://bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html) che abbiamo installato dal Bioconductor. Siamo fortunati a lavorare con Arabidopsis perché tra i circa 20 organismi supportati dal pacchetto abbiamo anche Ath. Questo pacchetto fornisce diverse possibilità e restituisce risultati sia tabulari che grafici. Inoltre permette di effettuare non solo arricchimenti ma anche classificazioni, sia per le categorie gene ontologies che per i pathway KEGG

# carichiamo la libreria
library(clusterProfiler)
library(org.At.tair.db)

# come primo comando andiamo ad esplorare le categorie alle quali appartengono i geni upregolati con un log2FC maggiore di 4
down <- res$log2FoldChange < 4
res[down,]
## log2 fold change (MLE): condition treated vs untreated 
## Wald test p-value: condition treated vs untreated 
## DataFrame with 23861 rows and 6 columns
##              baseMean log2FoldChange     lfcSE        stat       pvalue
##             <numeric>      <numeric> <numeric>   <numeric>    <numeric>
## AT1G01010  533.783591      1.1230292 0.4867571    2.307166 2.104559e-02
## AT1G01020  677.071094      0.6056971 0.2119271    2.858044 4.262606e-03
## AT1G01030  343.746971     -3.2819671 0.4509864   -7.277309 3.405457e-13
## AT1G01040  909.356189     -1.0993669 0.3179508   -3.457664 5.448807e-04
## AT1G01046    2.958823     -1.7600272 1.3278248   -1.325497 1.850065e-01
## ...               ...            ...       ...         ...          ...
## ATMG01370  168.494377     0.56781320 0.2883609  1.96910599   0.04894092
## ATMG01380   41.301320     0.06070519 0.6482423  0.09364583   0.92539050
## ATMG01390 4138.937078    -0.20154918 0.4538241 -0.44411299   0.65696090
## ATMG01400    2.616918     0.61534354 1.2322610  0.49936138   0.61752482
## ATMG01410   27.834388    -0.36224167 0.4516842 -0.80197988   0.42256460
##                   padj
##              <numeric>
## AT1G01010 4.039985e-02
## AT1G01020 9.764916e-03
## AT1G01030 3.323145e-12
## AT1G01040 1.506728e-03
## AT1G01046 2.627092e-01
## ...                ...
## ATMG01370   0.08431877
## ATMG01380   0.94502051
## ATMG01390   0.73263325
## ATMG01400           NA
## ATMG01410   0.51783617
down <- as.data.frame(res[res$log2FoldChange < -4,])
class(down)
## [1] "data.frame"
head(down)
##            baseMean log2FoldChange     lfcSE       stat       pvalue
## AT1G01250  34.87048      -4.027816 0.8262525  -4.874800 1.089187e-06
## AT1G01520   9.52135      -4.173945 1.0507842  -3.972219 7.120632e-05
## AT1G02205  86.68746      -4.972089 0.5546029  -8.965132 3.099096e-19
## AT1G02820 165.06660      -4.705437 0.3986524 -11.803356 3.750489e-32
## AT1G03310 641.44154      -5.226851 0.3371913 -15.501145 3.407656e-54
## AT1G04220  26.19707      -4.643998 0.8104425  -5.730200 1.003122e-08
##                   padj
## AT1G01250 4.773596e-06
## AT1G01520 2.326465e-04
## AT1G02205 5.348733e-18
## AT1G02820 1.538298e-30
## AT1G03310 4.469022e-52
## AT1G04220 5.819748e-08
dim(down)
## [1] 358   6
# come input di  groupGO usiamo la lista dei geni, che nel nostro caso corrisponde ai nomi di riga
# prima però dobbiamo convertire i TAIR in symbol e entrez
keytypes(org.At.tair.db)
##  [1] "ARACYC"       "ARACYCENZYME" "ENTREZID"     "ENZYME"      
##  [5] "EVIDENCE"     "EVIDENCEALL"  "GENENAME"     "GO"          
##  [9] "GOALL"        "ONTOLOGY"     "ONTOLOGYALL"  "PATH"        
## [13] "PMID"         "REFSEQ"       "SYMBOL"       "TAIR"
ids.down <- bitr(rownames(down), fromType="TAIR",
            toType=c("ENTREZID", "SYMBOL"),
            OrgDb="org.At.tair.db")
head(ids.down)
##        TAIR ENTREZID SYMBOL
## 1 AT1G01250   839322   <NA>
## 2 AT1G01520   839337   ASG4
## 3 AT1G02205   837602   CER1
## 4 AT1G02205   837602  CER22
## 5 AT1G02820   839304 AtLEA3
## 6 AT1G02820   839304   LEA3

CLASSIFICATION

ggo <- groupGO(gene     = ids.down$ENTREZID,
               OrgDb    =  org.At.tair.db,
               ont      = "BP",
               level    = 3,
               readable = T)
class(ggo)
## [1] "groupGOResult"
## attr(,"package")
## [1] "clusterProfiler"
# plot results
barplot(ggo, drop=TRUE, showCategory=20)

# transform to df to custom plot 
ggo.df <- as.data.frame(ggo)
head(ggo.df)
##                    ID                              Description Count
## GO:0019953 GO:0019953                      sexual reproduction     2
## GO:0019954 GO:0019954                     asexual reproduction     0
## GO:0022414 GO:0022414                     reproductive process    32
## GO:0032504 GO:0032504      multicellular organism reproduction     9
## GO:0032505 GO:0032505 reproduction of a single-celled organism     0
## GO:0061887 GO:0061887         reproduction of symbiont in host     0
##            GeneRatio
## GO:0019953     2/337
## GO:0019954     0/337
## GO:0022414    32/337
## GO:0032504     9/337
## GO:0032505     0/337
## GO:0061887     0/337
##                                                                                                                                                                                                                              geneID
## GO:0019953                                                                                                                                                                                                                  NA/DAA1
## GO:0019954                                                                                                                                                                                                                         
## GO:0022414 ATLOX3/ATLOX3/NA/DAA1/ATSUC5/ATSUC5/RAP2.7/RAP2.7/AtLEA4-2/AtLEA4-2/AtLEA4-2/DAD1/CYP78A6/CYP78A6/ATPAP16/ATPAP16/EMB1270/ATH1/NA/ATIPS1/ATIPS1/ATIPS1/ATIPS1/NA/ATNAS1/ATNAS1/AtLEA4-5/AtLEA4-5/WSD1/AtMYB24/AtMYB24/NA
## GO:0032504                                                                                                                                                                 ATLOX3/ATLOX3/NA/DAA1/AtLEA4-2/AtLEA4-2/AtLEA4-2/DAD1/NA
## GO:0032505                                                                                                                                                                                                                         
## GO:0061887
dim(ggo.df)
## [1] 585   5
str(ggo.df)
## 'data.frame':    585 obs. of  5 variables:
##  $ ID         : Factor w/ 585 levels "GO:0000075","GO:0000728",..: 160 161 174 207 208 449 406 49 100 101 ...
##  $ Description: Factor w/ 585 levels "abscission","acquisition of mycelium reproductive competence",..: 534 31 502 265 499 500 356 340 48 40 ...
##  $ Count      : int  2 0 32 9 0 0 14 139 31 176 ...
##  $ GeneRatio  : Factor w/ 40 levels "0/337","11/337",..: 14 1 18 38 1 1 7 6 17 8 ...
##  $ geneID     : chr  "NA/DAA1" "" "ATLOX3/ATLOX3/NA/DAA1/ATSUC5/ATSUC5/RAP2.7/RAP2.7/AtLEA4-2/AtLEA4-2/AtLEA4-2/DAD1/CYP78A6/CYP78A6/ATPAP16/ATPAP"| __truncated__ "ATLOX3/ATLOX3/NA/DAA1/AtLEA4-2/AtLEA4-2/AtLEA4-2/DAD1/NA" ...
# order levels to plot correctly
ggo.df$Description <- factor(ggo.df$Description,levels = ggo.df$Description[order(-ggo.df$Count)])

# order df and replace
ggo.df <- ggo.df[order(-ggo.df$Count),]

# plot first 50
ggplot(ggo.df[1:50,], aes(Description, Count)) + geom_point() +
  theme(axis.text.x=element_text(angle=90,hjust=1)) +
  ggtitle("GO most down regulated")

GO ENRICHMENT

ego <- enrichGO(gene          = ids.down$ENTREZID,
                OrgDb         = org.At.tair.db, 
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
                keyType       = 'ENTREZID',
                readable      = TRUE)
head(ego)
##                    ID                                          Description
## GO:0046148 GO:0046148                         pigment biosynthetic process
## GO:0009718 GO:0009718 anthocyanin-containing compound biosynthetic process
## GO:0042440 GO:0042440                            pigment metabolic process
## GO:0046283 GO:0046283    anthocyanin-containing compound metabolic process
## GO:0016144 GO:0016144                     S-glycoside biosynthetic process
## GO:0019758 GO:0019758                   glycosinolate biosynthetic process
##            GeneRatio   BgRatio       pvalue     p.adjust       qvalue
## GO:0046148    21/250 267/20467 1.585017e-11 1.111097e-08 9.676945e-09
## GO:0009718    11/250  60/20467 1.465387e-10 5.136180e-08 4.473286e-08
## GO:0042440    21/250 346/20467 1.882032e-09 4.397682e-07 3.830100e-07
## GO:0046283    11/250  85/20467 6.973849e-09 1.222167e-06 1.064430e-06
## GO:0016144    13/250 171/20467 1.920802e-07 1.923546e-05 1.675286e-05
## GO:0019758    13/250 171/20467 1.920802e-07 1.923546e-05 1.675286e-05
##                                                                                                                        geneID
## GO:0046148 NDF1/NA/AOR/ATMYB75/NA/SAUR6/ELIP/AT5MAT/GUN4/AtTT8/CCD4/ANS/ATH1/AtGLDP1/AtCHIL/ATGSTF12/DFR/UF3GT/PORA/LUT2/NdhN
## GO:0009718                                                    NA/ATMYB75/ELIP/AT5MAT/AtTT8/CCD4/ANS/AtCHIL/ATGSTF12/DFR/UF3GT
## GO:0042440 NDF1/NA/AOR/ATMYB75/NA/SAUR6/ELIP/AT5MAT/GUN4/AtTT8/CCD4/ANS/ATH1/AtGLDP1/AtCHIL/ATGSTF12/DFR/UF3GT/PORA/LUT2/NdhN
## GO:0046283                                                    NA/ATMYB75/ELIP/AT5MAT/AtTT8/CCD4/ANS/AtCHIL/ATGSTF12/DFR/UF3GT
## GO:0016144                               AtBE2/CYP79F2/BUS1/NA/FMO/ATLEUD1/ATGSTF11/NDF4/IPMI1/ATIPS1/ATMYB29/AtMYB76/ATMYB34
## GO:0019758                               AtBE2/CYP79F2/BUS1/NA/FMO/ATLEUD1/ATGSTF11/NDF4/IPMI1/ATIPS1/ATMYB29/AtMYB76/ATMYB34
##            Count
## GO:0046148    21
## GO:0009718    11
## GO:0042440    21
## GO:0046283    11
## GO:0016144    13
## GO:0019758    13
dotplot(ego)

#emapplot(ego)
cnetplot(ego, categorySize="pvalue")

#goplot(ego)
# provate se a voi vengono questi grafici

KEGG ENRICHMENT

kk <- enrichKEGG(gene         = rownames(down),
                 organism     = 'ath',
                 pvalueCutoff = 0.05)
head(kk)
##                ID                       Description GeneRatio BgRatio
## ath00941 ath00941            Flavonoid biosynthesis      4/66 22/4920
## ath00966 ath00966        Glucosinolate biosynthesis      4/66 23/4920
## ath00592 ath00592   alpha-Linolenic acid metabolism      5/66 42/4920
## ath01210 ath01210   2-Oxocarboxylic acid metabolism      5/66 74/4920
## ath00196 ath00196 Photosynthesis - antenna proteins      3/66 22/4920
##                pvalue    p.adjust      qvalue
## ath00941 0.0001802046 0.003888735 0.003411171
## ath00966 0.0002159520 0.003888735 0.003411171
## ath00592 0.0002160408 0.003888735 0.003411171
## ath01210 0.0029409165 0.031961771 0.028036641
## ath00196 0.0029594233 0.031961771 0.028036641
##                                                     geneID Count
## ath00941           AT4G22880/AT5G05270/AT5G42800/AT5G43935     4
## ath00966           AT1G16400/AT1G16410/AT2G43100/AT3G58990     4
## ath00592 AT1G17420/AT2G44810/AT3G03480/AT4G15440/AT5G38120     5
## ath01210 AT1G16400/AT1G16410/AT2G43100/AT3G02020/AT3G58990     5
## ath00196                     AT2G05070/AT2G05100/AT3G27690     3
dotplot(kk)

cnetplot(kk, categorySize="pvalue")

Come esercizio provate a ripetere le analisi di gene ontologies sulle restanti categorie (CC, MF) per i down-regolati e su tutte le categorie per gli up-regolati, considerando sempre un log2FC > 4